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ABSTRACT 

Radio-astronomical observations are increasingly corrupted by RF interference, and online 
detection and filtering algorithms are becoming essential. To facilitate the introduction of such 
techniques into radio astronomy, we formulate the astronomical problem in an array signal pro- 
cessing language, and give an introduction to some elementary algorithms from that field. We 
consider two topics in detail: interference detection by rank estimation of short-term covariance 
matrices, and spatial filtering by subspace estimation and projection. We discuss experimental 
data collected at the Westerbork radio telescope, and illustrate the effectiveness of the space- 
time detection and blanking process on the recovery of a 3C48 absorption line in the presence 
of GSM mobile telephony interference. 

Subject headings: methods: statistical; instrumentation: interferometers; methods: analytical 



1. INTRODUCTION 

Radio-astronomical observations are increasingly corrupted by RF interferers such as wireless com- 
munication and satellite navigation signals. Online detection and filtering algorithms are essential to reduce 
the effect of interference to an acceptable level. However, existing methods have a limited scope. Until 
now, the most widely implemented algorithm is a single-channel total power change detector, followed by 
a blanking of the correlator output. Friedman (1996) has implemented an improved power detector at the 
RATAN600, based on detection of change in the power. Weber et al. (1997) proposed the use of the 
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quantized correlation at all lags to test the presence of interference. Another detector based on wavelet de- 
composition has been proposed by Maslakovic et al. (1996). These are all single channel detectors which 
do not exploit the spatial properties of the interference. The only detector which considered combining mul- 
tiple telescopes for improved detection and blanking was proposed by Kasper Chute & Routledge (1982) 
for low frequency interferometry, where a robust data censoring method based on the temporal behavior of 
the cross spectrum was proposed. This requires a large number of estimated spectra (10 5 ) to obtain reliable 
robust estimates, and only two channels are used. Finally, adaptive filtering techniques have recently been 
considered by Barnbaum & Bradley (1998) who propose to excise interference from the Green-Bank radio 
telescope using a reference antenna and an LMS type algorithm. 

Our aim in this paper is to introduce modern array signal processing techniques to the context of radio 
astronomy, and to investigate the merits of multichannel detection and filtering algorithms at the Wester- 
bork Synthesis Radio Telescope (WSRT). By combining cross-correlation information of a large number of 
sensor pairs, we can increase the detection performance significantly, and also estimate the spatial signature 
of interferers. In essence, our approach is to compute (on-line) short-term spatial correlation matrices in 
narrow sub-bands, and then to compute the eigenvalue decomposition of each of these matrices (Leshem 
van der Veen & Deprettere 1999c). A rank estimate based on the eigenvalues allows to detect the number 
of interfering signals in each time-frequency slot, and the dominant eigenvectors give information on the 
"spatial signature" of the interferers. 

After detection, we can follow two directions. We can reduce the interference by rejecting corrupted 
time-frequency slots (blanking). This approach is suitable for time-slotted communication signals such as 
the European mobile telephony standard GSM, or the TDMA (time-division multiple access)-based mobile 
telephony standards IS -54/1 36 in the US. 

A more challenging approach is to also use the eigenvector information. Indeed, we can project out 
those dimensions in the spatial correlation matrices that correspond to the spatial signature vectors of the 
interference. Such spatial filtering techniques will greatly enhance the performance of observations with 
continuously -present interference. 

The effectiveness of the space-time detection and blanking process is demonstrated by applying the 
algorithms to data measured at the WSRT using an on-line 8-channel recording system. As will be shown 
in section 7, we were able to recover an absorption line of 3C48 which was completely masked by a super- 
imposed GSM interference, and could not be recovered by single channel techniques. 

The paper is written in a tutorial style, to appeal to both the radio astronomy and the signal processing 
communities. The structure of the paper is as follows. After posing the astronomical measurement equations 
in section 2, we reformulate the model in terms of array processing matrix language in section 3. We then 
introduce RF interference and describe its effect on the received data. In section 5 we discuss various 
detection algorithms. We compare the single and multichannel detectors, for the case of a narrow-band 
interferer with known spatial signature vector, and then present two multichannel detectors that do not 
assume this knowledge. We then move to spatial filtering techniques in section 6, where we formulate 
the basic ideas and describe a projections based approach. Finally, experimental results on multichannel 
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blanking are shown in section 7. 

2. ASTRONOMICAL MEASUREMENT EQUATIONS 

In this section we briefly describe a simplified mathematical model for the astronomical measurement 
process. Our discussion follows the introduction in Perley Schwab & Bridle (1989). The purpose of this is 
to connect to a matrix formulation commonly used in array signal processing, in the next section. 

The signals received from the celestial sphere may be considered as spatially incoherent wideband 
random noise. It is possibly polarized and perhaps contains spectral absorption or emission lines. Rather 
than considering the emitted electric field at a location on the celestial sphere, astronomers try to recover 
the intensity (or brightness) If (s) in the direction of unit-length vectors s, where / is a specific frequency. 
Let Ef(r) be the received celestial electric field at a location r on earth (see figure 1(a)). Assume that 
the telescopes are identical, and let A(s) denote the amplitude response of a telescope to a source in the 
direction s. The measured correlation of the electric fields between two sensors i and j with locations r-j and 
Tj is called a visibility and is (approximately) given by [Perley Schwab Sz Bridle (1989)] 

VfiTi^^EiEfMEfa)} = I A\ S )I f (s)e-^ sT ^y c dQ. 

J sky 

(E{ • } is the mathematical expectation operator, the superscript T denotes the transpose of a vector, and 
overbar denotes the complex conjugate.) Note that it is only dependent on the oriented distance r-j — rj 
between the two telescopes; this vector is called a baseline. 

For simplification, we may sometimes assume that the astronomical sky is a collection of d discrete 
point sources (maybe unresolved). This gives 

d 

/ /( S ) = J2lf(Sn)S{s-S n ), 
n=l 

where s ra is the coordinate of the n'th source, and thus 

d 

Vffarj) = A\s n )I f ( Sn )e-^ s ^- r ^ c . (1) 

n=l 

Up to this point we have worked in an arbitrary coordinate system. For earth rotation synthesis arrays, 
a coordinate system is often introduced as follows. We assume an array with telescopes that have a small 
field of view and that track a reference source direction so in the sky. Other locations in the field of view 
can be written as 

s = s + <x , S -L <T , 
(valid for small a) and a natural coordinate system is 

s = [0, 0, if, a = [£, m, 0] T . 
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Similarly, for a planar array, the receiver baselines can be parameterized as 

n - rj = X[u, v, w] T , A = j . 

The measurement equation in (u, v, w) coordinates thus becomes 

V M = ,-*"// A^ m)W , m) e-^U Mm . (2, 

The factor e -i 27 ™ is caused by the geometrical delay associated to the reference location, and can be com- 
pensated by introducing a slowly time-variant delay (see figure 1(6)). This synchronizes the center of the 
field-of-view and makes the reference source location appear as if it was at the north pole. After compensa- 
tion, we arrive at a measurement equation in (u, v) coordinates only, 



V f (u,v) = 




It has the form of a Fourier transformation. 

The function Vf(u, v) is sampled at various coordinates (u, v) by first of all taking all possible sensor 
pairs i,j or baselines r« — rj, and second by realizing that the sensor locations r-j, rj are actually time- 
varying since the earth rotates. Given a sufficient number of samples in the (u, v) domain, the relation can 
be inverted to obtain an image (the 'map'). 



3. ARRAY SIGNAL PROCESSING FORMULATION 

3.1. Obtaining the measurements 

We will now describe the situation from an array signal processing point of view. The signals received 
by the telescopes are amplified and downconverted to baseband. A time-varying delay for every telescope 
is also introduced, to compensate for the geometrical delay. 

Following traditional array signal processing practices, the signals at this point are called x.- L (t) rather 
than Ef(r), and are stacked in vectors 



x(t) 



xi(t) 
x p (t) 



where p is the number of telescopes. These are then processed by a correlation stage. 

It will be convenient to assume that x(i) is first split by a bank of narrow-band sub-band filters into a 
collection of frequency-components xj(t). The main output of the telescope hardware is then a sequence 
of empirical correlation matrices R/(t) of cross-correlations of xj(i), for a set of frequencies / € {/&} 
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covering a 10 MHz band or so, and for a set of times t € {i/J covering up to 12 hours. Each correlation 
matrix R/(t) is an estimate of the true covariance matrix R/(i), 

1 

R,(t) = E{ X/ (t) X/ (t) H } , R f (t) = -£ x /(* + nT )Mt + - (4) 

?1=0 

where the superscript H denotes a complex conjugate transpose, T is the sample period of x/(i) and N is the 
number of samples over which is averaged. This is drawn schematically in figure 2 (ignoring the detection 
stage for the moment). The matrices R/(£) are stored for off-line spectral analysis and imaging. 

Typically, the averaging period NT is in the order of 10-30 s, whereas each sub-band has a bandwidth 
in the order of 100 kHz or less. Due to the sub-band filtering, the original sampling rate of x(£) is reduced 
accordingly, resulting in T in the order of 10 /its. 

The connection of the correlation matrices R/(i) to the visibilities Vf(u, v) in section 2 is as follows. 
Each entry rjj(i) of the matrix R/(i) is a sample of this visibility function for a specific coordinate (u, v), 
corresponding to the baseline vector r-j(t) — rj(t) between telescopes i and j at time t: 

Ti(t) - Tj(t) = \[Uij(t), Vij(t), Wij(t)] 

Vf(uij(t),Vij(t)) = rij(t) . 

Note that we can obtain only a discrete set of (u, v) sample points. Indeed, the number of instantaneous 
independent baselines between p antennas is less than \p(p — 1). Also, using the earth rotation, the number 
of samples is given by the ratio of the observation time and the covariance averaging time (e.g., 12 h/30 
sec = 1440 samples). 

A few remarks on practical issues are in order. 

- Many telescope sites including WSRT follow actually a different scheme where the signals are first 
correlated at several lags and subsequently Fourier transformed. This leads to similar results. 

- The geometrical delay compensation is usually introduced only at the back end of the receiver. At this 
point, also a phase correction is needed to compensate for the factor e^ 27 ™^ W in the measurement equa- 
tion (2). This is referred to as fringe correction (Thompson et al. 1986). Since the earth rotates, Wij(t) 
is slowly time varying, with a rate of change in the order of 0-10 Hz depending on source declination 
and baseline length. 



3.2. Matrix formulation 

For the discrete source model, we can now formulate our measurement equations in terms of matrices. 
Let ro(ifc) be an arbitrary and time-varying reference point, typically at one of the elements of the array, and 
let us take the (u, v, w) coordinates of the other telescopes with respect to this reference, 

Ti(t) -r (t) = X[u i0 (t), v i0 (t), w i0 (t)} , i = !,-•• ,p. 
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Equation (1) can then be written slightly different as 

d 

V f {Ti{t),Tj(t)) = ^e-^ /s «( ri - r °)/ c ^ 2 (s n )/ / (s n )e J ' 27r/s " (r ^ ro)/c 



n=l 
d 



V f ( Uij (t), Vij (t)) = ^ e ~ jMuio{t)en+VM{t)mn)A ^n,rn n y 

n=l 

If(i n ,m n ) ■ A(£ n ,m n ) e J 2 ^Ko(t)^+%o(*)m„) _ 
In terms of correlation matrices, this equation can be written as 

R f (t) = A f (t)B f A H f (t) 

where 

A f (t) = [a ti /(4,mi), • • • ,& t j(e d ,m d )] 

and 

e -j(um(t)£+v m (t)m) 



(5) 



at,/(^,m) 



- j («p0 (t)^+WpO 



A(£ n , m Tl 



(6) 



B, = 



7/(4, mi) 







If{?d,m d ) 



The vector function a. t j(£,m) is called the array response vector in array signal processing. It describes 
the response of the telescope array to a source in the direction (£,m). As usual, the array response is 
frequency dependent. In this case, the response is also slowly time-varying due to the earth rotation. Note, 
very importantly, that the function as shown here is completely known, since the beam shape A(£, m) is 
calibrated and we know the locations of the telescopes very well. 

More realistically, the array response is less perfect. An important effect is that each telescope may 
have a different complex receiver gain, 7i(i), dependent on many angle-independent effects such as cable 
losses, amplifier gains, and (slowly) varying atmospheric conditions. If we take this into account, the model 
now becomes 



where 



R f (t) = T(t)A f (t)B f Ap)T(t) 

71 (*) 



r(t) = 



o 



o 



In future equations we will drop the dependence on /. 
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3.3. Additive noise 

In reality, most of the received signal consists of additive system noise. When this noise is zero mean, 
independent among the antennas (thus spatially white), and identically distributed, then it has a covariance 
matrix that is a multiple of the identity matrix, a 2 I, where a 2 is the noise power on a single antenna inside 
the subband which we consider. Usually the noise is assumed to be Gaussian. 

The resulting model of the received covariance matrix then becomes 

R(t) = r(t)A(t)BA H (t)r(t) H + a 2 I . (7) 

Note that this assumes that the noise is introduced after the varying receiver gains. This assumption is 
reasonable if the channels from the low-noise amplifier (LNA) outputs to the analog-to-digital converter 
(ADC) units are equal. Otherwise, it is still reasonable to assume that the noise is spatially white, i.e., 
the noise covariance matrix is diagonal. We can assume that it can be estimated using various calibration 
techniques; a simple diagonal scaling will then bring us back to the model (7). We further assumed that the 
quantization is fine, since a large dynamic range is needed to cope with strong interferers. 

The study of factorizations of the spatial covariance matrices such as shown above is the key to many 
array signal processing techniques. The knowledge of the specific structure of the array response vector (6) 
is of course already used in radio astronomy, as it enables the construction of the map using inverse Fourier 
techniques. The main point in this paper is to demonstrate that also interference adds a specific structure 
to the covariance matrices. This hopefully will allow its detection, provided our models are sufficiently 
accurate. 



4. RF INTERFERENCE 

RF interference (RFI) usually enters the antennas through the sidelobes of the main beam. It can be 
stronger or weaker than the system noise. An important property is that it has a certain spatial signature, or 
array response vector, which becomes explicit in the case of narrow-band signals. 

Examples of RFI present at WSRT are television broadcasts (Smilde station), geolocation satellites 
(GPS and Glonass), taxi dispatch systems, airplane communication and navigation signals, wireless mobile 
communication (GSM) and satellite communication signals (Iridium). Thus, interferers may be continuous 
or intermittent, narrow-band or wideband, and strong or weak. Some examples of actual interference are 
presented at the end of the section. 

Interference is usually not stationary over 10 seconds (let alone because of the time- varying fringe rate 
of 0-10 Hz), and hence it might be argued that it would average out from the long-term correlations. How- 
ever, the amount of nonstationarity is often insufficient to provide a good and reliable protection (Thompson 
1982), (Thompson et al. 1986). 
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4.1. Narrow-band interference model 



Suppose that we have a single interferer impinging on the telescope array. The interfering signal reaches 
the array with different delays r\ for each telescope. After demodulation to baseband, we have 



x i (t)=a l s(t-T i )e^ 2n ^ , 



i = l,- 



,P- 



Here, s(t) is the baseband signal, and cii represents the telescope gain in the direction of the interferer, 
including any possible attenuation of the channel. Unlike much of the array signal processing literature, 
the cij are likely to be different for each telescope since the interferer is typically in the near field. This 
implies that it impinges on each telescope at a different angle, whereas the response of the telescopes is not 
omni-directional . 

For narrow-band signals, time delays shorter than the inverse bandwidth amount to phase shifts of the 
baseband signal (Proakis 1995). This well-known property is fundamental to many phased array signal 
processing techniques. If the narrow-band assumption holds, then s(t — Tj) = s(t) and the model can be 
simplified. 

Note that we have already assumed before that the signals are sub-band filtered. Let W be the band- 
width of the sub-band filters. In WSRT, the largest baseline is 3000 m, corresponding to a maximal delay 
of 10 fis. Hence the narrow-band assumption holds if W <C 100 kHz Leshem van der Veen & Deprettere 
(1999c). Under this condition, we can stack the p telescope outputs from a particular sub-band filter in a 
vector x/(t) and write 



*/(*) 



-j2vr/Ti 



ape -j2nfr p 



s(t) =: as(t) 



As before, a is an array response vector. Unlike before, it is not a simple or known function of the direction 
of the interferer, since we are in the near field and the sidelobes of the array are not calibrated. The vector is 
also called the spatial signature of the interfering source. It is slowly time varying, and we write a = a(t). 



Similarly, with q interferers, 



q 



x/(t)=2^a J -(t) Sj (t) = A 8 (t)s(t), 



8(f) 



si(t) 



A s {t) = [ai(t),--- , 8i q (t)} . 



The subscript V is used to distinguish A s (t) from the array response matrix of the astronomical sources. 
The corresponding correlation matrix and its empirical estimate are 



M-\ 



R(t) = E{ X/ (t) X/ (t) H } = A a (t)R a (t)A?(t) , R(t) = ^ E x /(* + + mT ) . 



m=0 
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where R(t) is estimated by averaging over M samples. The q x (/-matrix R s (t) = E{s(t)s(t) H } depends on 
the correlation properties of the interfering signals. For independent interferers, it will be a diagonal matrix, 
with the q interfering powers on the diagonal. 

How well the estimate fits to R(i) depends on the stationarity of the scenario over the averaging in- 
terval MT, and is open to discussion. The power of television signals will be stationary over long periods 
(order tens of seconds or better). At the other extreme, communication signals such as used by the GSM 
mobile communication system are time slotted: time is partitioned into frames of about 5 ms and frames 
are partitioned into 8 slots. In this so-called time-division multiple access scheme (TDMA), each user can 
transmit only during its slot of 0.577 ms and then has to be silent for 7 times this period before transmitting 
again in the next frame. Thus, there is a short-term stationarity (over 0.577 ms), and a cyclostationarity with 
periods of about 5 ms. 

The stationarity of the columns of A s (t) depends on the stationarity of the location of the interferer, its 
distance, the fringe rate and the orientation of the telescopes. With multipath propagation, a mobile interferer 
only has to move about 30 cm to create a different a-vector, giving a stationarity in the order of 10-100 ms 
for a GSM user. Even for a fixed interferer such as a television station, the slow rotation of the telescopes 
as they track the sky will change the a-vector within a fraction of a second, either because of multipath 
fading or because the interferer moves through the highly variable sidelobe pattern. Another source of non- 
stationarity is the fringe correction introduced at the first IF stage to compensate for the geometrical delay. 
As the telescopes rotate, this introduces a time- varying phase, different to each telescope, with a rate in the 
range of — 10 Hz. 

The conclusion is that, for interference detection, R(t) is a useful estimate only over short averaging 
periods over which the interference is stationary, say MT in the order of milliseconds. Thus, M <C N, 
where NT x> 10 s, as in (4). 

4.2. Overall model: astronomical signals with interference and noise 

In summary, the model that we have derived is as follows: 

R(t) = T(t)A(t)BA iI (t)T(t) + A a (t)R a (t)A°(t) + a 2 I. 

R(t) is a p x p covariance matrix of which we have computed estimates at discrete times t. A : p x d is the 
array response matrix of the d discrete sources in the sky. Its columns are known functions of the (unknown) 
locations of the sources. It is a very wide matrix: d S> p, and assumed stationary over 10 s. B : d x d is 
a diagonal matrix (positive real) containing the brightness of each source, and assumed time-invariant over 
the complete observation. T are diagonal matrices (positive real) representing unknown and slowly varying 
antenna gains. 

A s : p x q is the array response matrix of the q interferers. It is likely to be unstructured. We will 
consider cases where q < p, so that A s is tall. H s : q x q is the interference correlation matrix. A s and R s 
are usually stationary only over very short time spans (order 10 ms). 
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a 2 I is the noise covariance matrix, assuming white independent and identically distributed noise for 
simplicity. The noise power a 2 is often rather well known. 

||ABA H ||, i.e., the observed power of the astronomical sources, is at least two orders of magnitudes 
smaller than a 2 , and for the purpose of detection, it can be ignored. In contrast, ||A S R S A^|| can be of 
comparable magnitude. 

4.3. Examples of interfering signals 

To demonstrate a few of its features, we present some measured observations of RFI. 1 

As mentioned before, interference may be continuous or intermittent. A prime example of continuously 
present interference are television broadcasts. Figure 3 shows a spectrogram centered at 780.75 MHz of the 
German television transmitter TV Lingen, located at about 80 km southeast of the WSRT. The two strong 
peaks in the spectrum are the two sound carrier waves. The received power of the TV station is much 
stronger than the WSRT system noise level, as can be seen from the fact that the baseband filter shape is 
barely visible. 

Figure 4 shows the GSM uplink band, which contains the communication of mobile handsets to the 
base stations. The short white dashes indicate the presence of (weak) interference. At least three channels 
can be seen at 902.4, 904.4, and 907.2 MHz, although there probably are more active channels at a lower 
power level. The TDMA time frame format of about 5 ms consisting 8 user slots of 0.577 ms can be 
recognized. Also visible is a weak continuous transmission at 902 MHz. This is likely to be an interference 
from the control building or an inter-modulation product. 

An example of the GSM downlink band (base station to mobiles) is shown in figure 5. Most of the 
signals are continuous in time, except for a few channels at e.g., 942.0, 942.8, and 949.8 MHz which are 
time slotted. 

Another interferer which one would like to remove from the observed data is the Iridium transmissions. 
Figure 6 shows a transmission of the Iridium satellite communication system at 1624 MHz (satellite-satellite 
communication and/or downlinks). It is clearly intermittent as well. 

Finally a wideband interferer is the GPS satellite navigation signal. This is a spread spectrum signal 
occupying a band of 10.23 MHz. Figure 7 shows a spectrogram of the GPS signal around 1575 MHz. 
One can clearly see the superposition of the narrow (1.023 MHz) commonly available C/A signal on the 
wideband (10.23 MHz) military P-code signal, resulting in the peak at the center frequencies. 



'The data has been collected using the NOEMI recording system described later in section 7. 
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5. INTERFERENCE DETECTION 



5.1. Introduction 



Ideally, the output of the correlation process produces clean estimates of A(t)BA(t) H , once every 10 s 
or so. In principle, we estimate it by 



As we have seen, these estimates are corrupted by interference and additive system noise, and unknown 
antenna gains. The objective of interference detection and rejection schemes is to improve the signal to 
interference and noise ratio (SINR) at the output of the integrators, i.e., at the 10 s level. Interference that 
is stationary at these time scales or longer can often be treated off-line. In this paper we consider online 
interference detection and excision schemes, assuming stationarity at millisecond time scales or less. 

Many interference detection schemes exist. They differ by the amount of knowledge that we can 
assume on the interfering signals, e.g., if we know the signal wave form, then the optimal detector has 
the form of a matched filter. Extensions are possible if the waveform is known up to a few parameters 
such as amplitude, phase or frequency. However, usually the signal is modulated by a message and hence 
effectively unknown. There are two classes of detection techniques: more or less deterministic methods that 
exploit known properties of the signals such as modulation type or certain periodicities, and those based 
on statistical models with unknown parameters, leading to Generalized Likelihood Ratio Tests (GLRT), a 
particular example of which is power detection. 

In principle, we can say that man-made interference is expected to be statistically different from the 
astronomical sources. Although this is a very attractive feature, it is not easy to use these properties for de- 
tection or excision, since the long averaging periods and the central limit theorem tend to jointly Gaussianize 
the interferers. For strong narrow-band interferers these methods are expected to give improved suppression 
at an increased computational expense (Leshem & van der Veen 1998). 

Another distinction between interferers and astronomical signals is their spatial signature vectors. As- 
tronomical signals enter through the main lobe of the telescopes and have a very structured (parametrically 
known) array response (viz. (6)), which is used for imaging. The interferers usually enter through the side- 
lobes and are in the near field, leading to unstructured a-vectors. Also, their location relative to the array 
is not correlated with the motion of the earth. It might even remain fixed relative to the array during the 
complete observation period (e.g., TV transmitters). Since the array tracks a fixed region in the sky which 
moves as the earth rotates, the directional vector of the interference is typically time varying. 

From all possibilities, we consider here two schemes: 

- Multichannel interference detection and excision. The interference is detected at short time scales (ms), 
and contaminated samples are removed from the averaging process in (8). This will work well if the 
interference is concentrated in frequency and time, as e.g., in the GSM system. 




N-l 



NT = 10 s. 



(8) 



n=0 



-12- 



- Spatial filtering. This more ambitious scheme is also suitable for continuously present interference such 
as TV stations. After detection, we estimate the spatial signature of the interferer and project out that 
dimension or otherwise subtract the signal coming from that direction. 

For the purpose of power detection schemes, it is sufficient to look at (short-term) correlation matrices 
based on measurement data in a window of length MT, with MT « 10 ms: 

^ M-l 

^ = mJ2 + mT ) x /(*fc + mT ) > = °> MT > 2MT ' ■ • • ■ 

m=0 

If an interferer is detected in this analysis window, it is discarded, otherwise the data is accepted and the 
correlation matrix is used in the formation of a clean estimate of Rj. 0s (i), as in figure 2. Obviously, many 
variations are possible, such as sliding window techniques, or discarding neighbors of contaminated samples 
as well (perhaps both in time and frequency). 

In this section we propose sub-band detection methods based on and analyze their performance. 
Spatial filtering is discussed in section 6. Throughout the section, we will drop the subscript k and write R 
and R for simplicity. 



5.2. Single channel spectral detector 

Detection theory is based on hypothesis testing. We test TLq\ there is no interference, versus Hi', there 
is at least one interferer in this band. The implementation of this test depends on the model that we pose for 
the interferer. We will first discuss some particularly simple cases which will allow analysis. 

Thus let us consider the single-channel case first. We assume that there is at most a single interferer, 
where the interfering signal is i.i.d. Gaussian noise with unknown power a 2 . The background noise is white 
Gaussian with known power a 2 . 

Without interferer, the observed data samples x, m = x(t m ) are complex normal (CM) distributed, with 
zero mean and variance a 2 . With an interferer, this distribution is still complex normal, but with variance 
a 2 + a 2 . Thus, we test the hypothesis 

Ho : x m ~ CAA(0, a 2 ) 

H\: x m ~ C7V r (0, a 2 + a 2 ) , m = 0, • • ■ , M - 1 . 
We assume that we have available M samples {x m }, collected in a vector x = [x\ , ■ ■ • , xm\- 

This is a rather standard problem in detection theory (see (Kay 1998) for an introduction). A Neyman- 
Pearson detector selects Hi if the likelihood ratio, 

_ P(x;fti) 

exceeds a threshold, where p(x; H) denotes the probability density function of x under the hypothesis H. 
It is known that this leads to an optimal probability of detection, given a certain probability of false alarm 
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(detecting an interferer when there is none). In our case, based on the model (9), the Neyman-Pearson 
detector simplifies to comparing the total received power to a threshold 7, deciding Hi if the test statistic 



^ M-l 

T(x) := ^ |^| 2 > 7. 

m=0 

Under the above assumptions we can obtain closed form expressions for the probability of false alarm and 
the probability of detection. For this, recall that the sum of squares of M real i.i.d. zero-mean unit- variance 
Gaussian random variables has a chi-square (x 2 ) distribution with M degrees of freedom. Since we have 
complex samples, T(x) is the sum-square of 2M real variables. Under Ho, these have a variance \, hence 
the probability of false alarm is given by 

P FA := P{T(x) > 7 ; H } = Q X | M (2 7 ) 

where Q X ? M (7) is the tail probability of a x 2 random variable with 2M degrees of freedom. It has a closed- 
form expression (cf. (Kay 1998)): 

M-i k 

fc=0 

Its inverse is known in terms of the inverse Gamma-function, and allows to select 7 to obtain a desired level 
of false alarm. Similarly, the probability of detection of an interference at this threshold 7 is given by 



P D := P{T(x) > 7 ; Hi} 

= ^E^=ikm| 2 >7; Hi} 



- a- 

M 



- p[ 2 T Is i 2 > 27 • H ^ (10) 

! xIm 1 1+INR 



Qy2^(i+inr) 



2 

where INR = % is the interference-to-noise ratio. 



5.3. Multichannel detector with known spatial signature 

A significant performance improvement is possible with a multichannel detector. To illustrate this, we 
assume again the simple case with at most a single narrow-band Gaussian interferer, with known spatial sig- 
nature vector a in white Gaussian noise. The source power of the interference is denoted by cr 2 ; to normalize 
the receiver gain we set ||a|| 2 := a H a = p, where p is the number of channels. Without interference, the 
data vectors x m are complex normal distributed with zero mean and co variance matrix cr 2 I. With a single 
interferer, the covariance matrix becomes R = E{x m x^} = cr 2 aa H + <r 2 I. Thus, 

H : x m ~£A/-(0,a 2 I) 

Hi: x m ~ CM(0, cr 2 aa H + cr 2 I) , m = 0, • • • , M - 1 . 
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The Neyman-Pearson detector based on the data matrix X = [xi , • • • , xm] considers the estimated data 
covariance matrix 

j M-l 
^ = ~M ^ X m X m 



M 

m=0 



and is given by (Kay 1998) 



r < x > : = ^sw | 7 ' 

This test is recognized as a matched spatial filter detector; essentially we compare the received energy in 
the direction a of the interferer to a 2 . If we define y rn to be the output of the matched beamformer in the 
direction of x m , 

_ 

Dm — I. n x m 

|a| 

then 

H : y m ~CN(0,a 2 ) 

Hi: y m ~ CM(0,pa 2 + a 2 ) , m = 0, • ■ ■ , M - 1 . 

and it is seen that taking the same threshold as in the single channel case will provide the same false alarm 
probability as before: 

P FA = P{T(X) > 7 ; H } = Q X | M (2 7 )- 
However, the probability of detection is now given by 

r° = rlT(X)>r, *,} - g^^JOLtf. 

Figure 8 presents the probabilities of detection as a function of interference to noise ratio for a single and for 
p = 14 channels. We have selected a threshold such that Pfa = 5%, which means that without interference, 
we will throw away 5% of the data. We can clearly see that the probability of detection is greatly improved 
by moving to the multichannel case. The improvement is equal to the array gain, 101og(p) = 11.5 dB. 



5.4. Single TDMA interferer with known spatial signature 

Let us now consider a TDMA signal: an interferer which is periodically active in a fraction j3 of the 
time (see figure 9). Here, < [3 < 1 is known as the duty cycle of the periodic signal. Assume that the 
interferer is present in the selected frequency band and that the duration of the slot in which the interferer 
is active is equal to aM samples x m , where we take a > 1. Let as before a\ denote the power of a single 
sample of the interferer when it is present. 
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Since the interfering slots need not be synchronized to the analysis window, a single interfering slot 
will give rise to two analysis windows in which the interferer is partially present, and possibly one or more 
analysis windows in which the interferer is present in all the samples. Since the interferer is time-slotted 
with duty cycle j3, there will also be windows that contain no interference. 

The corresponding probability density p(I) of having a certain average interference power / per sample 
in an arbitrary analysis window of length M can be computed in closed form, as 

a + 1 



p(I) 



(1 



a 



-P) 5(1) , 



1 2 



a 



P, 



a 



a 







< I < I„ 



It is plotted in figure 9, where the vertical arrows indicate the unit impulse function S( ■ ). For example, for 
an interferer of strength a 2 s per sample when it is on, the maximal average interference power per sample is 
obviously of, when all samples are contaminated. The probability of this is (a — l)/a (3. Power densities 
less than of occur with a uniform distribution for analysis windows that are only partly corrupted, at the 
edges of the interference slot. 



We can define 



- the average interference power per sample before detection: 




- the average interference power per sample after detection and blanking: 

Ires = J I(l-P D (I))p(I)dI, 

- the fraction of number of samples kept after detection and blanking: 

n res = J (l-P D (I))p(I)dI. 

Figure 10 shows the dependence of the residual INR as a function of M (the number of samples in an 
analysis block), for an interferer of length L = 64 sub-band samples, a duty cycle = 1/8, and a false alarm 
rate of 5%. Obviously, very weak interference is not detected, and in that case we throw away 5% of the 
data due to the false alarm rate. High interference powers are easily detected, and almost all contaminated 
analysis windows will be detected and blanked. Only the tails of an interfering slot might be missed, so that 
there is still some interference remaining after detection. The worst case occurs for interference that is not 
strong enough to be detected all the time, but not weak enough to be harmless. 



Several other interesting facts can be seen in these figures. The most important is the large performance 
gain in the multichannel approach, as compared to a single channel. As seen in figure 8, the effect of using 
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an array is to shift the graphs of probability of detection to the left by the array gain, e.g., for the 14-channel 
detector the graph is shifted by 11.5 dB. Hence, we require 11.5 dB less interference power in order to detect 
it. However, the effective gain is given by the vertical distance between the graphs: this shows the amount of 
interference suppression for a given interference power. In figure 10 the suppression can be approximately 
21 dB larger than that of the single antenna case. 

A second interesting phenomenon is the fact that the interference suppression is almost the same for 
a large range of analysis windows M. Thus, we would take this window rather small, so that the residual 
number of samples is larger. This effect is mainly due to the fact that the case of partial blocks with weaker 
power is less frequent as the analysis block becomes shorter. Further study of this model appeared in 
(Leshem k, van der Veen 1999). 

5.5. Eigenvalue analysis 

So far, we have looked at the detection problem from a rather idealistic viewpoint: at most 1 interferer, 
and a known spatial signature. The reason was that for this case, we could derive optimal detectors with 
closed-form expressions for the performance. We will now discuss an extension to more practical situations. 

Our goal is the detection of the presence of an interferer from observed correlation data. As a start, let 
us first consider the covariance matrix due to q interferers and no noise, 

R = A S R S A" 

where R has size p x p, A s has size p x q and R s has size q x q. For a low number of interferers q, this 
brings us to familiar grounds in array signal processing, as it admits analysis by subspace-based techniques. 
We give a brief introduction here; see (Krim Sz Viberg 1996) for an overview and references. 

If q < p, then the rank of R is q since A s has only q columns. Thus, we can estimate the number of 
narrow-band interferers from a rank analysis. This is also seen from an eigenvalue analysis: let 

R = UAU H 

be an eigenvalue decomposition of R, where the p x p matrix U is unitary (UU H = I, U H U = I) and 
contains the eigenvectors, and the p x p diagonal matrix A contains the corresponding eigenvalues in non 
increasing order (Ai > A2 > • • • > A p > 0). Since the rank is q, there are only q nonzero eigenvalues. We 
can collect these in a q x q diagonal matrix A s , and the corresponding eigenvectors map x q matrix U s , 
so that 

R = U S A S U*. (11) 

The remaining p — q eigenvectors from U can be collected in a matrix U n , and they are orthogonal to XJ S 
since U = [U s U n ] is unitary. The subspace spanned by the columns of XJ S is called the signal subspace, 
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the orthogonal complement spanned by the columns of U n is known as the noise subspace (although this is 
a misnomer). Thus, in the noise-free case, 



R = UAU H = [U s U„] 















In the presence of white noise, 



R = A s R s A* + a 2 I E 



U H 

^ s 

U H 



(I p denotes apxp identity matrix.) In this case, R is full rank: its rank is always p. However, we can still 
detect the number of interferers by looking at the eigenvalues of R. Indeed, the eigenvalue decomposition 
is derived as (expressed in terms of the previous decomposition (11) and using the fact that U = [U s U n ] 
is unitary: U S U^ + U„U* = I p ) 



R = A s R s A^ + cj 2 I p 

= XJ S A S U^ + a 2 (U S U^ + U n U£) 

= U s (A s + a%)U^ + U n (a%_ g ; 

= [U s U n ] 

=: UAU H 



A s + a 2 I q 








cr 2 lp-q 



(12) 



hence R has p — q eigenvalues equal to a 2 , and q that are larger than a 2 . Thus, we can detect the number of 
interferers q by comparing the eigenvalues of H to a threshold defined by a 2 . 

A physical interpretation of the eigenvalue decomposition can be as follows. The eigenvectors give 
an orthogonal set of "directions" (spatial signatures) 2 present in the covariance matrix, sorted in decreasing 
order of dominance. The eigenvalues give the power of the signal coming from the corresponding directions, 
or the power of the output of a beamformer matched to that direction. Indeed, let the i'th eigenvector be u;, 
then this output power will be 



ufRu, 



A,; 



The first eigenvector, ui , is always pointing in the direction from which most energy is coming. The second 
one, U2, points in a direction orthogonal to ui from which most of the remaining energy is coming, etcetera. 

If there is no interference and only noise, then there is no dominant direction, and all eigenvalues are 
equal to the noise power. If there is a single interferer with power a 2 and spatial signature a, normalized to 
||a|| 2 = p, then the covariance matrix is R = cr 2 aa H + a 2 I. It follows from the previous that there is only 



2 Here, direction is not to be interpreted as the physical direction-of-incidence of the interferer, but rather the abstract direction 
of a unit-norm vector in the vector space C p . Due to multipath, unequal gains and fringe corrections, the physical direction-of- 
incidence might not be identifiable from the spatial signature a. 
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one eigenvalue larger than a 2 . The corresponding eigenvector is ui = a-pyy, and is in the direction of a. 
The power coming from that direction is 

Ai = UjRuj = pa 2 + a 2 . 

Since there is only one interferer, the power coming from any other direction orthogonal to ui is a 2 , the 
noise power. Note the connection with the test statistic of the previous section, where we assumed that a is 
known. Since ui = a-p^, 

a H Ra _ uf Rui _ 

H — H — ^1 • 

a a u x ui 

Thus, the test statistic of the previous section reduces to testing the dominant eigenvalue of R, and knowl- 
edge of a is in fact not needed. 

With more than one interferer, this generalizes. Suppose there are two interferers with powers o\ and 
CT2, and spatial signatures ai and a2. If the spatial signatures are orthogonal, afa2 = 0, then ui will be in the 
direction of the strongest interferer, number 1 say, and Ai will be the corresponding power, Ai = pa 2 + a 2 . 
Similarly, A2 = pa 2 + a 2 . 

In general, the spatial signatures are not orthogonal to each other. In that case, ui will point into the 
direction that is common to both ai and a2, and 112 will point in the remaining direction orthogonal to ui. 
The power Ai coming from direction ui will be larger than before because it combines power from both 
interferers, whereas A2 will be smaller. 

The covariance matrix eigenvalue structure can be nicely illustrated on data collected at the WSRT. We 
selected a narrow band slice (52 kHz) of a GSM uplink data file, around 900 MHz. In this subband we have 
two interfering signals: a continuous narrow band CW signal which leaked in from a local oscillator, and 
a weak GSM signal. From this data we computed a sequence of short term cross spectral matrices R,0- 5ms 
based on 0.5 ms averages. Figure 11 shows the time evolution of the eigenvalues of these matrices. The 
largest eigenvalue is due to the CW signal and is always present. The GSM interference is intermittent: at 
time intervals where it is present the number of large eigenvalues increases to two. The remaining eigenval- 
ues are at the noise floor, a 2 . The small step in the noise floor after 0.2 s is due to a periodically switched 
calibration noise source at the input of the telescope front ends. 

The eigenvalue decomposition (12) shows more than just the number of interferers. Indeed, the columns 
ofU s span the same subspace as the columns of A s . This is clear in the noise-free case (11), but the 
decomposition (12) shows that the eigenvectors contained in U s and U n respectively are the same as in the 
noise-free case. Thus, 

span(U s ) = span(A s ) , U^A S = . (13) 

Given a correlation matrix R estimated from the data, we compute its eigenvalue decomposition. From this 
we can detect the rank q from the number of eigenvalues larger than a 2 , and we can determine XJ S and hence 
the subspace spanned by the columns of A s . Although we cannot directly identify each individual column 
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of A s , its subspace estimate can nonetheless be used to filter out the interference — such spatial filtering 
algorithms are discussed in section 6. Note that it is crucial that the noise is spatially white. For colored 
noise, an extension (whitening) is possible but we have to know the coloring. 



5.6. Multichannel detector with unknown spatial signature 

In case we only have an estimate R based on a finite amount of samples M and the spatial signature 
vectors of the interference are unknown, there are no optimal results. The eigenvalue analysis suggested that 
we should compare the eigenvalues to a threshold defined by a 2 : without interference, all eigenvalues are 
asymptotically equal to a 2 . We will discuss two detectors, one for the case where a 2 is known, and another 
one for which it is unknown. 

If the noise power a 2 is known, we can apply the (generalized) likelihood ratio test (GLRT), which 
leads to a method due to Box (1949) for testing the null hypothesis that a R = I (no interference). The 
GLRT leads to a test statistic given by 

v ~ 

T(X) := -MplogJJ-i (14) 

i=l 

where A, is the 2-th eigenvalue of R, and we detect an interferer if T(X) > 7. This basically tests if all 
eigenvalues are equal to a 2 , with a certain confidence. In the no-interference case, one can show that 

nx) ~ 4+i )(P -2) 

This allows to select the value of 7 to achieve a desired false alarm rate. 

If also the noise power is unknown, we propose to use the Minimum Description Length (MDL) detec- 
tor [Wax Sz Kailath (1985)]. In this case, rather than setting a threshold based on the asymptotic distribution 
of the LRT, we try to find the correct model order which minimizes the description length of the data. The 
MDL rank estimator is given by 

q = argminMDL(n) (15) 

n 

where 

( nf=n+i ^» ) p 1 

MDL(n) = -(p -n)M log— x — + -n(2p - n + 1) log M 

p—n Si=n+1 

and an interference is detected if q 7^ 0. The first term basically tests if the geometric mean of the smallest 
p — n eigenvalues is equal to the arithmetic mean, which is only true if these eigenvalues are equal to each 
other. (The second term is a correction that grows with the number of unknown parameters to be estimated). 
Note also that the arithmetic mean of the small eigenvalues is an estimate of the noise variance, so in the 
case of testing whether n = or not the first term in the MDL reduces to a sample estimate T(x) of (14). 
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This rank detector is simple to implement since it is independent of the varying SINR in the system. A 
disadvantage is that the false alarm rate is not known and not fixed. 

Finally a simple option which can be used to limit the false alarm rate is to collect a number of process- 
ing blocks, sort them according to the value of the statistic T(x), defined in (14) and throw away a given 
percentage with the highest score. This is conceptually simple but needs more memory available. 

Experimental results on multichannel blanking are presented in section 7. 

6. SPATIAL FILTERING 

Let us now assume that we have obtained a covariance matrix R, which contains the rather weak 
covariance matrix of the astronomical sources (visibilities) Rj,, plus white noise. Suppose also that there is 
an interferer with power a 2 : 

R = + ofaa 11 + a 2 I . 

In the previous section, we considered schemes to detect the interferer from the eigenvalues of R, a short- 
term estimate of R. After detection, we proposed to discard R from a longer-term average if it is found to 
be contaminated, but what if the interferer is present all the time? In that case, it is more suitable to try to 
suppress its contribution a 2 &a H . This leads to spatial filtering techniques. 

6.1. Projecting out the interferer 

An elementary form of spatial filtering is to null all energy with spatial signature a. To this end, we can 
introduce the p x p projection matrix 

P^ = I-a(a H a)' 1 a H . 

Pa is a projection because Pa = Pa- ^ * s a l so easily seen that P^a = 0: this direction is projected 
out. If we denote by R the filtered covariance matrix, we obtain 

R := PgRPf = ParR.PaT + a 2 P^ . (16) 

Thus, the interference is removed by the projection. At the same time, the visibility matrix is modified by 
the projections, and the noise is not white anymore, since one dimension is missing. The imaging stage has 
to be aware of this, which is the topic of (Leshem & van der Veen 1999b). 

In general, a is not known. However, note that we do not need a, but only a projection matrix to project 
it out. Recall from equation (12) the eigenvalue decomposition of R, and in particular the matrix containing 
an orthonormal basis of the "noise subspace" U n , which is the orthogonal complement of a, with p — 1 
columns. According to (13), U^a = 0. It is now straightforward to show that 

P^ = U n U* (17) 
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Indeed, since U"U ra = I p -i, 

p-Lp-L _ TJ TJ H TJ TJ H — TJ TJ H — p-L 

and 

P^a = U n U^a = 0. 

Thus, we can compute the required projection matrix directly from the eigenvalue decomposition of R. 

Expression (17) can immediately be generalized to the more general case of q < p interferers and 
unknown a-vectors. Indeed, in this case, the projection onto the complement of the A s -matrix of the inter- 
ference is given by 

Pi s = I - A^A,)-^ = UnK 

Note that we do not have to know A s : the relevant noise subspace is estimated from the eigenvalue decom- 
position of R. This hinges upon the fact that the noise covariance is white (in general: known), and the 
visibility matrix Rt, is insignificant at these time scales (otherwise, it might disturb the eigenvalue decom- 
position). 

As an alternative to (16), we can define another filtered covariance matrix 

R := U*RU n = U^R,U n + a 2 l^ q , (18) 

where we have used U ra _L A s and U^U n = I p - q . In this case, R has size (p — q) x (p — q). Al- 
though smaller, this matrix contains the same information as P^RP^. Besides the dimension reduction, an 
advantage of this scheme is that the noise term stays white. 



6.2. Keeping track of projections 

Since the projections alter the visibility data in R„, it is essential, for the purpose of imaging, to 
store the linear operation represented by the projections. At the same time, it might be necessary to adapt 
the projection several times per second, since the a-vectors of interferers are time-varying. Hence, in the 
construction of the 10 s correlation average from short-term projected correlation matrices, we also have to 
construct the effective linear operation. 

Let Rfc denote the short-term correlation matrix, where k = 0, 1, ■ • ■ , — 1 is the time index, and N is 
the number of short-term matrices used in the long-term average. Denote for generality the linear operation 
representing the projection by where = (U ra )fc(U n )" in the first filtering scheme (equation (16)), 
and Lfc = (U n )^ in the second (equation (18)). Consider now the short-term filtered averages, 



R fe :— L fc R fc Lf — LfcRfLf + <r 2 L fc Lf , k = 0, 1, • • • , N — 1 . 
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By simply averaging these, the long-term average will be 

^ JV-1 j N-l 



N 



k=0 k=0 

The appear here at both sides of R&. To move them to one side, we make use of the general expression 

vec(ABC) = (C T ® A)vec(B) 
where <g> denotes a Kronecker product, and vec(-) the column- wise stacking of a matrix into a vector, 



A<8>B 



anB ai 2 B 

02lB 022B 



A = [ai a 2 • • • ] 



vec(A) := 



ai 
a 2 



In this case, we obtain 



vec(R 10s ) 



= jf E [(_ L fc ® L fe )vec(R fe )] 

= [iE L t® L fc] vec(R„) + o- 2 [jf L fc ® L fc] vec(Ip) 
= Cvec(R„) + cr 2 Cvec(Ip) 



where 



1 

C — L k$ L /.- 



fc=0 

and the overbar denotes complex conjugation. C is the effective linear mapping of entries of H v to entries 
of R 10s . For the imaging step, we have to know how R 10s depends on H v . Thus, we have to construct and 
store C along with R 10s . Unfortunately, it is a rather large matrix: p 2 x p 2 in the first filtering scheme, and 
(p — q) 2 x p 2 in the second. Another problem for imaging might be that the noise contribution on R 10s is 
no longer white, but determined by C. Two possible remedies are 

- Assume that the a-vectors were sufficiently variable over the time interval. In that case, C is likely to be 
of full rank and thus invertible, and we can construct 



C" i vec(R 10s ) = vec(R^) + a\ec{I p ) . 



By unstacking the result, we recover our interference-free model H v + a 2 1. However, the inversion of 
C might be a formidable, and numerically dubious, task. 
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- If we take = (U n )£ as in (18), then the noise contribution on each R& is white. We can average the 
Rfc if they have the same dimension, i.e., p — q where the number of interferers q is constant over the 
interval. In that case, 

1 N ^ 

n=0 

so that the noise contribution on R 10s is white. Note that no inversion is necessary. 

If we do not invert C then the observed visibilities V(uij, Vij) in the matrix H v are modified by some 
(known) linear combination. This has implications for the synthesis imaging step. The usual construction 
of an image using inverse Fourier transformation (based on (3)) now gives rise to a point-source image 
convolved with a space-varying point spread function ("dirty beam"). Since the point spread function is 
known at every location in the image, it is still possible to correct for it using an extension of the usual 
CLEAN deconvolution algorithm. Details are in (Leshem & van der Veen 1999b). 

Since C is a factor p 2 larger than R 10s , it might in fact be more efficient to store the sequence of spatial 
filters Lfc. This is the case if Lfc is to be updated at time scales of 10 s/p 2 = 50 ms or longer. 



6.3. Other spatial filtering possibilities 

Without going into too much detail, we mention a few other possibilities for spatial filtering and inter- 
ference cancellation. Suppose there is a single interferer, 

R = R„ + o- 2 aa H + a 2 I . 

- Subtraction. With an estimate of a and a 2 , we can try to subtract it from the covariance data: 

R = R - ct 2 M h . (19) 

Without other knowledge, the best estimate of a is the dominant eigenvector, m, of R, and likewise the 
best estimate of a 2 is Ai — a 2 . Since both of these are derived from R, it turns out to be not too different 
from the projection scheme. Indeed, if we look at 

(I — auiuf )R(I — craiUi) = R — uiu"Ai(2q — a 2 ) 

we can make it equal to (19) by selecting a such that Ai(2a — a 2 ) = a 2 . The projection scheme had 
a = 1. 

Our point here is that also subtraction can be represented by a two-sided linear operation on the corre- 
lation matrix. Consequently, the visibility matrix R„ is altered, and hence the corrections mentioned in 
section 6.2 are in order. 
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- Subtraction of a reference signal. If we have a reference antenna that receives a 'clean' copy of the 
interfering signal, then we might try to subtract this reference signal from the telescope signals. There 
are many adaptive schemes for doing so, e.g., the LMS algorithm Haykin (1995). The general scheme is 
as illustrated in figure 12. In this figure, the a- vector of the interferer is found by cross-correlating with 
the reference antenna. We also estimate its power. After correcting for the noise power on the reference 
antenna, we can reconstruct and subtract &s(t). 

This scheme is rather similar to the original projection approach where we reduce the dimension to the 
noise subspace, viz. equation (18). The main difference is that, now, we reduce the dimension fromp + 1 
antennas back to p antennas, so there is no loss of dimensions from the astronomy point of view. It 
appears that this only has advantages if the reference antenna has a better INR than the telescopes. Also, 
we need as many reference antennas as there are interferers. 

As with the projection technique, all of these forms of spatial filtering modify the observed visibilities 
in the matrix Rt, by a known linear combination, with implications for the synthesis imaging step (Leshem 
& van der Veen 1999b). 

7. MULTICHANNEL BLANKING: EXPERIMENTAL RESULTS 

To test the blanking and filtering algorithms, we have attached the WSRT antennas to a multi-channel 
data recorder that can collect a few seconds of data at 20 MHz rate and store it on CDROM. This enabled 
us to record a variety of actual interference and process it off-line. In this section, we demonstrate the 
performance of the blanking algorithm by adding GSM observations to "clean" galactic 3C48 data, in a 
variety of scalings. The results are quite good, as it is possible to recover a 3C48 absorption line which was 
completely masked by the GSM interference. 

7.1. Experimental setup 

The data recorder has been acquired in the context of the STW NOEMI project, a cooperation between 
Delft University of Technology and ASTRON/NFRA. It basically consists of an industrial PC with four 
PCI.212 sampling boards. Each board contains two ADCs, and the boards are synchronized so that in 
total eight telescope channels can be sampled simultaneously. The ADCs have a resolution of 12 bit with 
sampling rates of 20 MHz down to 0.313 MHz in steps of a factor of 2. After collecting a batch of data, it 
can be copied into system memory (384 MB), previewed and stored onto CDROM. 3 

Fig. 13 shows an overview of the WSRT system to indicate the point where the NOEMI data recorder 
was connected. The WSRT is an East- West linear array of fourteen telescope dishes, mostly spaced at 
144 m. Each dish is equipped with a front-end receiver that can be tuned to several frequency bands. 



3 We would like to thank G. Schoonderbeek for programming the data acquisition software. 
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Both polarizations (X and Y) are received. The resulting 14 x 2 channels are amplified, filtered, down- 
converted to an intermediate frequency (IF) range around 100 MHz, and transported to the main building 
via coaxial cables. Here, the signals are fed to the equalizer unit which compensates for the frequency 
dependent attenuation in the ground cables. The equalizer unit has outputs for the broadband continuum 
system (DCB, 8 bands of 10 MHz) and for the spectral line system (DLB, 10 MHz). In the equalizer unit 
and in the DCB/DLB IF systems are mixers, amplifiers and filter units which take care of the baseband 
conversion and filtering. At baseband the signals are digitized to 2-bit resolution, a correction is applied 
for the geometrical delay differences between the telescopes, and the signals are correlated (in pairs) in the 
DZB/DCB correlators. The NOEMI recorder is connected at the output of the DLB spectral line IF system. 
Of the 14 x 2 available telescope channels, a selection of eight are connected to the NOEMI ADC samplers. 
The geometrical delay compensations and fringe rate corrections were not included in the measurements. 4 

The WSRT system contains also calibration noise sources, which are switched on for a 1.25s period 
every 10 seconds. For regular WSRT observations these noise sources are used for system noise and gain 
calibration purposes. In some of the observed NOEMI data sets these noise sources are clearly visible as a 
5-15% power step. 

Two important tests have been applied to the recording system. The synchronization of the channels 
was checked by applying a common wideband signal and was found to be in order. The cross-talk between 
the channels was measured by applying a signal to only one of the channels and looking at the leakage into 
the other channels. The power insulation between two channels on the same PCI board is found to be 51 dB 
(0.28% in voltage), and at least 90 dB (0.0032% in voltage) across boards. This is sufficient for spectral line 
work and for RFI mitigation tests. 



7.2. Clean 3C48 absorption data 

To compare our off-line frequency domain correlation process based on recorded data to the online 
Westerbork correlator (the DZB) we have made an interference-free observation of the galactic HI absorption 
of 3C48, a spectral line at 1420 MHz. Figure 14 shows the estimate of the power spectral density of the 
received signal based on the largest eigenvalue of the covariance matrix. 

The coherency (correlation coefficient) of signals x.- L and xj at the output of telescopes i and j is defined 

as 

_ mom = . m m 

Since all telescopes are tracking the same source s, we have that x\ = ajS + rii where n,i is the noise at 



4 Since these fringe rates are in the order of 0-10 Hz, this has no consequences for the detection of interference based on 
short-term correlation matrices, with typical integration periods in the order of milliseconds. 
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telescope i. With uncorrelated noise of power E(|nj| 2 ) = a 2 , and a source power of a 2 , it follows that 

Pij{f) = ai " j a 2 J) + * 2 (f) {l + 3Y 

Thus, the theoretical value of the coherency is constant over all nonzero baselines, and can be estimated 
based on the parameters of 3C48 and knowledge of the receiver gains and noises. These theoretically 
expected (asymptotic) values can then be compared to the computed coherencies of the recorded observation 
using (20), and can also be compared to the coherency measured with the DZB hardware. 

Figure 15(a) shows the magnitude of the coherency function for all nonzero baselines as based on a 
NOEMI recording of a few seconds. The coherency is around 5%, and the spectral absorption at 1420.4 MHz 
shows up as a dip. We verified that the absorption line is statistically significant. For comparison we include 
the same spectral line as processed by the WSRT DZB correlator in figure 15(6). The values of the coherency 
are in good agreement (differences are due to differences in processing bandwidths, observation times and 
instrumental settings). 

To verify the phase behavior of the coherency we have computed the unwrapped phase as a function 
of frequency. Note that the geometrical delay compensation and fringe corrections were not included in the 
recording. Due to the narrowband processing, the delay offset Tij of one channel with respect to another 
shows up as a frequency-dependent phase shift e _j,27r ^ T ^ (the fringe), which will be the phase of Pij(f). 
Here, depends on the location of 3C48 and the baseline vector r-j — rj between antenna i and j, and 
is known. Figure 16 compares the observed phase differences (averaged over all identical baselines) to the 
computed phase, as a function of frequency and baseline length. It is seen that the correspondence is very 
good. Note that for the shorter baselines we have more realizations so that their correspondence is better. 



7.3. 3C48 absorption line with GSM interference 

At this point we are ready to demonstrate the performance of the sub-band detection and blanking 
method as presented in section 5. To this end, we have superimposed on the 3C48 data (at 1420 MHz) 
another measurement file containing GSM interference (at 905 MHz), with the same bandwidth and for 
various amplitude scalings of each file. Although a bit artificial, the good linearity of the WSRT system 
implies that had a GSM signal been transmitted with a carrier frequency of 1420 MHz, then the measured 
data would be the superposition of the two signals plus system noise. The overlay allows us to verify 
the blanking performance for various mixtures of signal-to-interference power, since the clean data is now 
available as a reference and also the theoretical coherency is well known. 

As described in section 5, the detection of an interferer in a specific time-frequency cell is based on the 
eigenvalues of the corresponding correlation matrix of the resulting mixture. In this scheme, if one or more 
eigenvalues are above a threshold, then an interferer is detected and that data block is omitted. However, to 
avoid the selection of the threshold based on a desired false alarm rate, we have chosen to simply throw away 
the worst 30 percent of the data according to the value of the detector. We have computed the coherency 
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of the clean, the contaminated and the blanked signals. Figure 17 shows the coherency functions over all 
baselines for a particular mixture of signals and interference: scaling the GSM data file by 0.1 and the 
clean 3C48 data file by 0.9. It is seen that (a) the clean 3C48 spectrum shows the absorption line, which is 
(b) completely masked when GSM interference is added. After blanking, (d) the absorption line is almost 
perfectly recovered. For comparison we also included (c) the results of blanking based on single channel 
power detection from channel 2 only, without the sub-band decomposition. The failure of this common way 
of single channel detection is clearly seen. The reason is that the GSM signal was rather weak, so that for 
single-channel wideband processing the probability of detection was quite low, even for a false alarm rate 
ofupto30%. 

To show the effect of interference power we have repeated the experiment with the GSM data set 
weighted by a factor 0.5. The stronger GSM interferer is now more easily detected and the resulting spec- 
trum after blanking is better as seen in figure 18. 

8. Conclusions 

In this paper, we considered various aspects of multichannel interference suppression for radio-astronomy. 
It was shown that by sub-band processing we have access to the many narrow-band techniques available in 
array signal processing. We have demonstrated the benefits of multichannel spatio-spectral blanking, both 
theoretically and on measured data. The results are very pleasing. We have also discussed spatial filter- 
ing techniques and demonstrated how they can be incorporated into the radio-astronomical measurement 
equation. 
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Fig. 1. — (a) The emitted electrical field from the celestial sphere is received by a rotating telescope array; 
(b) geometrical delay compensation 
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Fig. 2. — The astronomical correlation process 
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Fig. 3. — Television broadcast 
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Fig. 5.— GSM downlink 
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Fig. 7. — GPS transmission, showing the civil code (BW= 1.023 MHZ) superimposed on the wideband 
military code (BW= 10.23 MHZ). f c = 1575 MHz. 
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P D vs. INR 
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Fig. 8.— P D vs. INR, for M = 10, 30, 64, and false alarm rate P FA = 5% 
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Fig. 9. — (a) Interferer with slot length L = aM samples, power a 2 s per on-sample, and duty cycle (3. (b) 
Corresponding probability density of interference power in a single analysis window. 
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Fig. 10. — (a) Effective residual INR after blanking versus effective INR at the input; (b) fraction of remain- 
ing samples after blanking 




Fig. 11. — Eigenstructure as a function of time 
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12. — Estimation of a using a reference antenna 
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13. — Overview of main WSRT systems with the NOEMI data recorders 
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Fig. 14. — 3C48: largest eigenvalue of the covariance matrix 




Fig. 15. — 3C48 coherency function, magnitude (for all baselines), (a) NOEMI recording and off-line 
processing, (b) online WSRT processing by the DZB 
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Fig. 16. — 3C48 averaged coherency phase function vs. frequency, various baselines. 
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Fig. 17. — Magnitude of the coherency functions of 3C48 mixed with GSM interference. GSM data scaled 
by 0.1. (a) clean 3C48 data, (b) 3C48 mixed with GSM, (c) after single channel detection/blanking, (d) 
after multichannel subband detection/blanking 
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3C48 Coherency: No interference With GSM interference 
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Fig. 18. — Magnitude of the coherency functions of 3C48 mixed with GSM interference. GSM data scaled 
by 0.5. (a) clean 3C48 data, (b) 3C48 mixed with GSM, (c) after single channel detection/blanking, (d) 
after multichannel subband detection/blanking 



